// --*- C++ -*------x-----------------------------------------------------------
//
// Author:          Enrico Negri, Alessandro Albiero, Silvio Tosatto
//
// Project name:    BiocompUP Align
//
// Date:            01/2008
//
// Description:     A program for checking the manipulation of profiles.
//
// -----------------x-----------------------------------------------------------

#include <Alignment.h>
#include <AlignmentBase.h>
#include <GetArg.h>
#include <IoTools.h>
#include <Profile.h>
#include <iostream>


using namespace Biopool;


/// Show command line options and help text.
void
sShowHelp()
{
	cout << "\nPROFILE CHECKING"
	     << "\nA program for checking the manipulation of profiles.\n"
	     << "\nOptions:"
	     << "\n * [--in <name>] \t Name of input BLAST mode 6 file"
	     << "\n   [-d <double>] \t Min. master vs all seq. id. filter threshold (default = 0.00)"
	     << "\n   [-D <double>] \t Min. all vs all seq. id. filter threshold (default = 0.00)"
	     << "\n   [-u <double>] \t Max. master vs all seq. id. filter threshold (default = 1.00)"
	     << "\n   [-U <double>] \t Max. all vs all seq. id. filter threshold (default = 1.00)"
	     << "\n   [--fasta]     \t Use FASTA instead of BLAST M6"
	     << "\n   [--gap]       \t Allow gaps in target sequence"
	     << "\n   [--verbose]   \t Verbose mode"
	     << "\n" << endl;
}


int
main(int argc, char **argv)
{
	double downs, downa, ups, upa;
	string inputFileName;


	// --------------------------------------------------
	// 0. Treat options
	// --------------------------------------------------

	if (getArg("h", argc, argv))
	{
		sShowHelp();
		return 1;
	}

	getArg("-in", inputFileName, argc, argv, "!");
	getArg("d", downs, argc, argv, 999.9);
	getArg("D", downa, argc, argv, 999.9);
	getArg("u", ups, argc, argv, 999.9);
	getArg("U", upa, argc, argv, 999.9);
	bool fasta = getArg("-fasta", argc, argv);
	bool gap = getArg("-gap", argc, argv);
	bool verbose = getArg("-verbose", argc, argv);

	string examplesPath = getenv("VICTOR_ROOT");

	if (examplesPath.length() < 3)
		cout << "Warning: environment variable VICTOR_ROOT is not set." << endl;
	string path = getenv("VICTOR_ROOT");
	examplesPath = path + "examples/";
	
	Alignment ali;
	if (inputFileName != "!")
	{
		inputFileName = examplesPath + inputFileName;
		ifstream inputFile(inputFileName.c_str());
		if (!inputFile)
			ERROR("Error opening input BLAST mode 6 file.", exception);
		if (fasta)
			ali.loadFasta(inputFile);
		else
			ali.loadBlastMode6(inputFile);
	}
	else
		ERROR("checkprofile needs input BLAST mode 6 file.", exception);


	// --------------------------------------------------
	// 1. Edit alignment
	// --------------------------------------------------

	if (downs < 999.9)
	{
		cout << "Setting simple lower seq. id. threshold at " << downs << endl;
		ali.RemoveLowerSimple(downs);
	}
	else
		if (downa < 999.9)
		{
			cout << "Setting full lower seq. id. threshold at " << downa << endl;
			ali.RemoveLowerAll(downa);
		}

	if (ups < 999.9)
	{
		cout << "Setting simple upper seq. id. threshold at " << ups << endl;
		ali.RemoveUpperSimple(ups);
	}
	else
		if (upa < 999.9)
		{
			cout << "Setting full upper seq. id. threshold at " << upa << endl;
			ali.RemoveUpperAll(upa);
		}

	Profile profile;
	if (gap)
		profile.setAllowGaps(true);
	profile.setProfile(ali);


	// --------------------------------------------------
	// 2. Print profile data
	// --------------------------------------------------

	fillLine(cout);
	cout << "\n     #     Gaps     Max AA     Max Freq     Sum Freq\n" << endl;
	for (unsigned int i = 0; i < profile.getSequenceLength(); i++)
	{
		cout << setw(6)  << i
		     << setw(9)  << profile.getNumGap(i)
		     << setw(11) << profile.getAminoMaxFrequency(i)
		     << setw(13) << setprecision(3) << profile.getFreqMaxAminoFrequency(i);

		double tmp = 0.00;
		for (AminoAcidCode j = ALA; j <= TYR; j++)
			tmp += profile.getAminoFrequencyFromCode(j, i);
		cout << setw(13) << setprecision(3) << tmp << endl;
	}
	cout << endl;
	fillLine(cout);

	if (verbose)
	{
		cout << "\nMaster sequence:\n"        << profile.getSeq()
		     << "\n\nConsensus sequence:\n"   << profile.getConsensus()
		     << "\n\nLenght of profile    = " << profile.getSequenceLength()
		     << "\nSequences in profile = "   << profile.getNumSequences()
		     << "\n" << endl;
		fillLine(cout);
	}

	return 0;
}
